Data from isolated Langerhans islets with Chromium Next GEM Single Cell 5p RNA library v2. Aligned with STARsolo.
library(Seurat)
library(Matrix)
fpath <- '/home/vqf/proyectos/TFGs/Ruben/NoteBooks/mat/alignments/'
fsuffix <- 'Solo.out/Gene/filtered'
smpls <- c('158WT', '238WT', '244WT', '161KO', '235KO', '243KO')
types <- c('WT', 'WT', 'WT', 'KO', 'KO', 'KO')
tmpobjs <- c()
for (i in (1:length(smpls))){
smpl <- smpls[i]
tpe <- types[i]
mtx <- paste(fpath, smpl, '/', smpl, fsuffix, '/matrix.mtx', sep = '')
cll <- paste(fpath, smpl, '/', smpl, fsuffix, '/barcodes.tsv', sep = '')
ftr <- paste(fpath, smpl, '/', smpl, fsuffix, '/features.translated.tsv', sep = '')
interm <- ReadMtx(mtx = mtx, cells = cll, features = ftr)
langer <- CreateSeuratObject(interm, project = 'Langerhans')
langer$orig.ident <- smpl
langer@meta.data$genotype <- tpe
tmpobjs <- c(tmpobjs, langer)
rm(interm)
}
langer <- merge(x = tmpobjs[[1]], y = tmpobjs[2:length(tmpobjs)], add.cell.ids=smpls)
str(langer)
Formal class 'Seurat' [package "SeuratObject"] with 13 slots
..@ assays :List of 1
.. ..$ RNA:Formal class 'Assay5' [package "SeuratObject"] with 8 slots
.. .. .. ..@ layers :List of 6
.. .. .. .. ..$ counts.1:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. .. .. .. .. .. ..@ i : int [1:12625864] 132 133 195 280 396 407 542 549 625 648 ...
.. .. .. .. .. .. ..@ p : int [1:9136] 0 1306 5962 7093 9605 10537 11822 13335 14320 16832 ...
.. .. .. .. .. .. ..@ Dim : int [1:2] 56980 9135
.. .. .. .. .. .. ..@ Dimnames:List of 2
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. ..@ x : num [1:12625864] 1 1 1 1 1 1 1 1 2 3 ...
.. .. .. .. .. .. ..@ factors : list()
.. .. .. .. ..$ counts.2:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. .. .. .. .. .. ..@ i : int [1:12764382] 450 620 1088 1270 1533 1629 1656 2337 2492 2730 ...
.. .. .. .. .. .. ..@ p : int [1:9347] 0 209 970 1626 2549 3625 5719 6952 8045 8578 ...
.. .. .. .. .. .. ..@ Dim : int [1:2] 56980 9346
.. .. .. .. .. .. ..@ Dimnames:List of 2
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. ..@ x : num [1:12764382] 1 1 1 1 1 1 1 1 1 1 ...
.. .. .. .. .. .. ..@ factors : list()
.. .. .. .. ..$ counts.3:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. .. .. .. .. .. ..@ i : int [1:8332109] 149 239 251 455 495 556 558 625 628 712 ...
.. .. .. .. .. .. ..@ p : int [1:9225] 0 1187 2246 3194 3866 4500 5387 5856 6938 7669 ...
.. .. .. .. .. .. ..@ Dim : int [1:2] 56980 9224
.. .. .. .. .. .. ..@ Dimnames:List of 2
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. ..@ x : num [1:8332109] 1 1 1 2 1 1 1 1 1 5 ...
.. .. .. .. .. .. ..@ factors : list()
.. .. .. .. ..$ counts.4:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. .. .. .. .. .. ..@ i : int [1:21186415] 125 185 239 251 271 319 481 490 503 556 ...
.. .. .. .. .. .. ..@ p : int [1:12284] 0 1235 2980 5220 6477 8583 9381 11054 12627 12875 ...
.. .. .. .. .. .. ..@ Dim : int [1:2] 56980 12283
.. .. .. .. .. .. ..@ Dimnames:List of 2
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. ..@ x : num [1:21186415] 1 1 1 4 1 1 2 2 1 3 ...
.. .. .. .. .. .. ..@ factors : list()
.. .. .. .. ..$ counts.5:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. .. .. .. .. .. ..@ i : int [1:22457424] 87 149 150 172 233 239 251 271 274 278 ...
.. .. .. .. .. .. ..@ p : int [1:11602] 0 1566 2561 4913 5904 8185 9762 11143 12669 19210 ...
.. .. .. .. .. .. ..@ Dim : int [1:2] 56980 11601
.. .. .. .. .. .. ..@ Dimnames:List of 2
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. ..@ x : num [1:22457424] 1 1 1 1 2 1 3 7 1 1 ...
.. .. .. .. .. .. ..@ factors : list()
.. .. .. .. ..$ counts.6:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. .. .. .. .. .. ..@ i : int [1:22697606] 22 149 206 228 229 251 301 322 380 427 ...
.. .. .. .. .. .. ..@ p : int [1:10124] 0 3062 5843 7576 9831 15634 21951 24022 27154 29339 ...
.. .. .. .. .. .. ..@ Dim : int [1:2] 56980 10123
.. .. .. .. .. .. ..@ Dimnames:List of 2
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. ..@ x : num [1:22697606] 1 1 1 1 1 2 1 3 2 2 ...
.. .. .. .. .. .. ..@ factors : list()
.. .. .. ..@ cells :Formal class 'LogMap' [package "SeuratObject"] with 1 slot
.. .. .. .. .. ..@ .Data: logi [1:61712, 1:6] TRUE TRUE TRUE TRUE TRUE TRUE ...
.. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. .. .. .. ..$ : chr [1:61712] "158WT_AAACCTGAGTGGTCCC" "158WT_AAACCTGCAATGCCAT" "158WT_AAACCTGCACCCATGG" "158WT_AAACCTGCAGATCTGT" ...
.. .. .. .. .. .. .. ..$ : chr [1:6] "counts.1" "counts.2" "counts.3" "counts.4" ...
.. .. .. .. .. ..$ dim : int [1:2] 61712 6
.. .. .. .. .. ..$ dimnames:List of 2
.. .. .. .. .. .. ..$ : chr [1:61712] "158WT_AAACCTGAGTGGTCCC" "158WT_AAACCTGCAATGCCAT" "158WT_AAACCTGCACCCATGG" "158WT_AAACCTGCAGATCTGT" ...
.. .. .. .. .. .. ..$ : chr [1:6] "counts.1" "counts.2" "counts.3" "counts.4" ...
.. .. .. ..@ features :Formal class 'LogMap' [package "SeuratObject"] with 1 slot
.. .. .. .. .. ..@ .Data: logi [1:56980, 1:6] TRUE TRUE TRUE TRUE TRUE TRUE ...
.. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. .. .. .. ..$ : chr [1:56980] "Gm37671" "Gm19087" "Gm8941" "Gm38212" ...
.. .. .. .. .. .. .. ..$ : chr [1:6] "counts.1" "counts.2" "counts.3" "counts.4" ...
.. .. .. .. .. ..$ dim : int [1:2] 56980 6
.. .. .. .. .. ..$ dimnames:List of 2
.. .. .. .. .. .. ..$ : chr [1:56980] "Gm37671" "Gm19087" "Gm8941" "Gm38212" ...
.. .. .. .. .. .. ..$ : chr [1:6] "counts.1" "counts.2" "counts.3" "counts.4" ...
.. .. .. ..@ default : int 6
.. .. .. ..@ assay.orig: chr(0)
.. .. .. ..@ meta.data :'data.frame': 56980 obs. of 0 variables
.. .. .. ..@ misc : list()
.. .. .. ..@ key : chr "rna_"
..@ meta.data :'data.frame': 61712 obs. of 4 variables:
.. ..$ orig.ident : chr [1:61712] "158WT" "158WT" "158WT" "158WT" ...
.. ..$ nCount_RNA : num [1:61712] 3918 24759 2859 7331 1551 ...
.. ..$ nFeature_RNA: int [1:61712] 1306 4656 1131 2512 932 1285 1513 985 2512 912 ...
.. ..$ genotype : chr [1:61712] "WT" "WT" "WT" "WT" ...
..@ active.assay: chr "RNA"
..@ active.ident: Factor w/ 1 level "Langerhans": 1 1 1 1 1 1 1 1 1 1 ...
.. ..- attr(*, "names")= chr [1:61712] "158WT_AAACCTGAGTGGTCCC" "158WT_AAACCTGCAATGCCAT" "158WT_AAACCTGCACCCATGG" "158WT_AAACCTGCAGATCTGT" ...
..@ graphs : list()
..@ neighbors : list()
..@ reductions : list()
..@ images : list()
..@ project.name: chr "SeuratProject"
..@ misc : list()
..@ version :Classes 'package_version', 'numeric_version' hidden list of 1
.. ..$ : int [1:3] 5 1 0
..@ commands : list()
..@ tools : list()
rm(tmpobjs)
counts1 <- GetAssayData(object = langer, layer = "counts.1")
counts1[200:210, 1:3]
11 x 3 sparse Matrix of class "dgCMatrix"
158WT_AAACCTGAGTGGTCCC 158WT_AAACCTGCAATGCCAT 158WT_AAACCTGCACCCATGG
Epha4 . . .
Mir6352 . . .
Colgalt2 . . .
Gm28386 . . .
Gm28385 . . .
Gm26263 . . .
Gm28387 . . .
Poglut2 . . .
Btf3-ps14 . . .
Catspere2 . . .
Gm37313 . . .
Now all objects are in langer. Starting sample is in
langer@meta.data$orig.ident and genotype is in
langer@meta.data$genotype.
table(langer@meta.data$genotype)
KO WT
34007 27705
langer <- PercentageFeatureSet(langer, pattern = "^mt-", col.name = "percent.mt")
table(langer@meta.data$genotype)
KO WT
34007 27705
Visualize channels for claim
VlnPlot(langer, features = c("nFeature_RNA"), group.by = "orig.ident", pt.size = 0, y.max = 7500)
VlnPlot(langer, features = c("nCount_RNA"), group.by = "orig.ident", pt.size = 0, y.max = 7500)
VlnPlot(langer, features = c("percent.mt"), group.by = "orig.ident", pt.size = 0)
VlnPlot(langer, features = c("nFeature_RNA", "nCount_RNA","percent.mt"), group.by = "orig.ident", ncol = 3, pt.size = 0)
VlnPlot(langer, features = c("nFeature_RNA", "nCount_RNA","percent.mt"), group.by = "genotype", ncol = 3, pt.size = 0)
reads_per_cell <- langer[['nFeature_RNA']]$nFeature_RNA
MIN_READS_PER_CELL <- 300
MAX_READS_PER_CELL <- 6000
MAX_MT <- 5
plot(sort(reads_per_cell), xlab='cell', log='y', main='Reads per cell (ordered)')
abline(h=MIN_READS_PER_CELL, col='magenta')
abline(h=MAX_READS_PER_CELL, col='blue')
table(langer@meta.data$genotype)
KO WT
34007 27705
langer <- subset(langer, percent.mt<MAX_MT)
langer.highexpr <- subset(langer, nFeature_RNA>MAX_READS_PER_CELL)
langer <- subset(langer, nFeature_RNA>MIN_READS_PER_CELL)
langer <- subset(langer, nFeature_RNA<MAX_READS_PER_CELL)
VlnPlot(langer, features = c("nFeature_RNA", "nCount_RNA","percent.mt"), group.by = "orig.ident", ncol = 3, pt.size = 0)
VlnPlot(langer, features = c("nFeature_RNA", "nCount_RNA","percent.mt"), group.by = "genotype", ncol = 3, pt.size = 0)
langer <- SCTransform(langer, vars.to.regress = "percent.mt", verbose = FALSE)
table(langer@meta.data$genotype)
KO WT
31588 26893
langer <- RunPCA(langer)
DimHeatmap(langer, dims = 1:10, cells = 50, balanced = TRUE)
ElbowPlot(langer)
langer <- RunTSNE(langer, dims = 1:20)
langer <- RunUMAP(langer, dims = 1:20)
Using method 'umap'
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
TSNEPlot(langer, split.by="genotype")
FeaturePlot(langer, features = c("Gcg", "Ins1", "Ins2", "Jchain", "Sst", "Ppy", "Ghrl", "Cd68", "Cd4", "Pecam1", "Tagln", "Iapp", "Col3a1"), cols = c("grey", "red"), reduction = "tsne", combine = F, split.by = "orig.ident")
[[1]]
[[2]]
[[3]]
[[4]]
[[5]]
[[6]]
[[7]]
[[8]]
[[9]]
[[10]]
[[11]]
[[12]]
[[13]]
[[14]]
[[15]]
[[16]]
[[17]]
[[18]]
[[19]]
[[20]]
[[21]]
[[22]]
[[23]]
[[24]]
[[25]]
[[26]]
[[27]]
[[28]]
[[29]]
[[30]]
[[31]]
[[32]]
[[33]]
[[34]]
[[35]]
[[36]]
[[37]]
[[38]]
[[39]]
[[40]]
[[41]]
[[42]]
[[43]]
[[44]]
[[45]]
[[46]]
[[47]]
[[48]]
[[49]]
[[50]]
[[51]]
[[52]]
[[53]]
[[54]]
[[55]]
[[56]]
[[57]]
[[58]]
[[59]]
[[60]]
[[61]]
[[62]]
[[63]]
[[64]]
[[65]]
[[66]]
[[67]]
[[68]]
[[69]]
[[70]]
[[71]]
[[72]]
[[73]]
[[74]]
[[75]]
[[76]]
[[77]]
[[78]]
VlnPlot(langer, features = c("Gcg", "Ins1", "Ins2", "Gapdh", "Ghrl"), group.by = "genotype", pt.size = 0)
langer <- FindNeighbors(langer, reduction = 'pca', dims = 1:20)
langer <- FindClusters(langer, resolution = 0.5)
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 58481
Number of edges: 2037766
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.9553
Number of communities: 28
Elapsed time: 12 seconds
DimPlot(langer, reduction = 'tsne', split.by = 'genotype', label = T, repel = T)
Cluster 17 probably contains beta and delta cells ###Get the genes used to classify the clusters
options(future.globals.maxSize= 8288490189)
langer <- PrepSCTFindMarkers(langer, assay = "SCT", verbose = FALSE)
all.markers <- FindAllMarkers(langer)
View(all.markers)
saveRDS(langer, "/home/vqf/proyectos/TFGs/Ruben/NoteBooks/langer.rds")
saveRDS(all.markers, "/home/vqf/proyectos/TFGs/Ruben/NoteBooks/allmarkers.rds")